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Abstract 

Composites  with  extremal  or  unusucd  thermal  expansion  coefficients  are  designed 
using  a  three-phase  topology  optimization  method.  The  composites  are  made  of  two 
different  material  phases  and  a  void  phase.  The  topology  optimization  method  consists 
in  finding  the  distribution  of  material  phases  that  optimizes  an  objective  function  (e.g., 
thermoelastic  properties)  subject  to  certain  constraints,  such  as  elastic  symmetry  or 
volume  fractions  of  the  constituent  phases,  within  a  periodic  base  ceU.  The  effective 
properties  of  the  material  structures  are  found  using  the  numericd  homogenization 
method  based  on  a  finite-element  discretization  of  the  base  cell.  The  optimization 
problem  is  solved  using  sequential  linear  programming. 

To  benchmark  the  design  method  we  first  consider  two-phase  designs.  Our  optimal 
two-phase  microstructures  are  in  fine  agreement  with  rigorous  bounds  and  the  so-called 
Vigdergauz  microstructures  that  realize  the  bounds.  For  three  phases,  the  optimal 
microstructures  are  also  compared  with  new  rigorous  bounds  and  again  it  is  shown 
that  the  method  yields  designed  materials  with  thermoelastic  properties  that  are  close 
to  the  bounds. 

The  three-phase  design  method  is  illustrated  by  designing  materials  having  max¬ 
imum  directional  thermal  expansion  (thermal  actuators),  zero  isotropic  thermal  ex¬ 
pansion,  and  negative  isotropic  thermal  expansion.  It  is  shown  that  materials  with 
effective  negative  thermal  expansion  coefficients  can  be  obtained  by  mixing  two  phases 
with  positive  thermal  expansion  coefficients  and  void. 


•Partly  on  leave  from  Department  of  Solid  Mechanics,  Technical  University  of  Denmark,  DK-2800  Lyngby, 
Denmark  (current  and  permanent  address). 
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1  Introduction 


In  this  paper,  we  use  a  topology  optimization  procedure  to  determine  the  distribution  of  three 
phases  (two  different  bulk  material  phases  and  a  void  phase)  in  order  to  design  composites 
with  extremal  or  unusual  thermal  expansion  behavior.  Three  phases  are  used  (as  opposed 
to  two  phases)  since  one  can  achieve  effective  properties  of  the  composite  beyond  those  of 
the  individual  components  [25].  Microstructural  variation  is  limited  to  one  length  scale  in  a 
unit  cell  as  this  is  most  easily  manufacturable. 

Materials  with  extreme  or  unusual  thermal  expansion  behavior  are  of  interest  from  both  a 
technological  and  fundamental  standpoint.  Of  particular  practical  interest  are  materials  with 
zero  thermal  expansion,  maximum  thermal  expansion  or  force,  and  negative  (i.e.,  minimum) 
thermal  expansion.  Heretofore,  however,  a  systematic  procedure  to  design  materials  with 
exotic  thermal-expansion  behavior  has  been  lacking. 

Materials  with  zero  thermal  expansion  coefficients  are  needed  for  use  in  structures  subject 
to  temperature  changes  such  as  civil  engineering  and  space  structures  as  well  as  piping 
systems.  Examples  are  bridges,  where  temperature  changes  between  day  and  night,  and 
summer  and  winter,  cause  big  structural  changes,  and  space  applications,  where  temperature 
differences  between  sunny  and  shady  sides  of  a  structure  are  extreme.  Such  temperature 
differences  can  cause  distortion  of  space  antennas  and  "rapid”  temperature  changes  due  to 
orbit  of  the  Hubble  space  telescope  are  known  to  cause  thermal  distortions  of  its  solar  arrays 
in  turn  causing  the  arrays  and  thereby  the  telescope  to  vibrate  (Collins  and  Richter  [8]). 
Materials  with  maximum  unidirectional  thermal  displacement  or  force  can  be  employed  as 
“thermal”  actuators.  Materials  with  negative  thermal  expansion  coefficients  can  be  used  to 
overcome  positive  thermal  expansion  of  other  materials  or,  among  other  applications,  be  used 
for  thermally  operated  fasteners.  A  fastener  made  of  a  negative  thermal  expansion  coefficient 
material,  upon  heating,  can  be  inserted  easily  into  a  hole  because  of  the  volume  contraction. 
When  cooled  down,  it  will  expand,  fitting  tightly  into  the  hole  and  upon  heating  can  be 
easily  removed.  Finally,  the  design  of  material  with  specific  thermal  expansion  coefficients 
is  important,  to  be  able  to  eliminate  thermal  mismatch  between  parts  in  structures  subject 
to  heat  changes. 

A  negative  thermal  expansion  material  has  the  counterintuitive  property  of  contracting 
upon  heating.  There  are  a  number  of  existing  materials  with  negative  thermal  expansion  co¬ 
efficients.  Glasses  in  the  titania-silica  family  have  isotropic  negative  expansion  coefficients  at 
room  temperature  [40].  Examples  of  materials  that  have  negative  thermal  expansion  at  very 
low  temperatures  (<  100/f)  are  silicon  and  germanium  [22]  as  well  as  Bi2.2Sri_sC aCu20x 
superconductor  single  crystals  [48].  Examples  of  materials  with  directional  negative  thermal 
expansion  coefficients  at  room  temperature  are  Kevlar,  carbon  fibers,  plastically  deformed 
(anisotropic)  Invar  (Fe-Ni  alloys)  [18]  and  certain  molecular  crystalline  networks  [3].  The 
negative  expansion  mechanism  of  these  molecular-level  networks  is  based  on  un-twisting  of 
helical  chains.  Currently  there  is  no  way  to  manufacture  these  materials  in  extended  form, 
but  this  an  active  area  of  research. 

An  interesting  question  [3]  is  whether  there  is  a  mechanistic  relationship  between  negative 
thermal  expansion  and  negative  Poisson’s  ratio?  A  material  with  negative  Poisson’s  ratio 
expands  laterally  when  pulled  axially  and  can  be  manufactured  by  processing  of  open-walled 
foam  structures  described  by  Lakes  [24].  We  will  show  that  isotropic  materials  with  effective 
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negative  thermal  expansion  coefficients  exist  with  positive  values  of  the  Poisson’s  ratio  and 
that  they  can  be  obtained  by  mixing  two  phases  with  positive  thermal  expansion  coefficients 
and  void. 

Several  researchers  have  addressed  the  problem  of  designing  materials  composites  with 
specific  directional  thermal  expansion  properties.  Autio,  Laitinen  and  Pramila  [2]  designed 
laminates  with  specific  elastic  and  thermal  expansion  coefficients  by  varying  layering  thick¬ 
nesses  and  directions.  Wetherhold  and  Wang  [47]  tailored  the  thermal  deformation  of  beams, 
and  Parton  and  Kudryavtsev  [31]  discussed  the  design  of  one  dimensional  beams  with  nega¬ 
tive  thermal  expansion.  Rodriques  and  Fernandes  [34]  designed  thermally  loaded  structures 
with  optimal  stiffness. 

It  was  shown  by  Levin  [27]  and  later  by  Rosen  and  Hashin  [35],  that  there  is  a  sim¬ 
ple  relationship  between  the  effective  thermal  expansion  coefficients  and  the  effective  elastic 
moduli  of  two-phase  materials.  In  other  words,  designing  two-phase  composites  with  extreme 
thermal  expansion  coefficients  corresponds  to  designing  two-phase  composites  with  extreme 
bulk  moduli.  The  problem  of  finding  the  structures  that  extremize  the  effective  elastic 
properties  of  two-phase  media  has  a  long  history  beginning  with  the  composite-sphere  as¬ 
semblages  of  Hashin  and  Shtrikman  [17]  for  the  bulk  modulus  problem.  Certain  hierarchical 
laminates  were  shown  to  realize  the  Hashin-Shtrikman  bounds  on  both  the  bulk  and  shear 
moduli  of  isotropic  two-phase  composites  [10].  More  recently,  Milton  and  Cherkaev  [28]  have 
found  multi-length  scale  materials  possessing  elastic  properties  ranging  over  the  entire  range 
compatible  with  thermodynamics.  Vigdergauz  [45,  46]  and  Grabovsky  and  Kohn  [12,  13] 
have  studied  single-inclusion  microstructures  of  extreme  rigidity.  Sigmund  [41,  42,  43]  has 
designed  material  structures  with  specific  elastic  properties  (including  isotropic  negative 
Poisson’s  ratio  material),  where  the  microstructure  is  restricted  to  one  length  scale. 

For  three-phase  materials,  one-to-one  relationships  between  the  thermal  expansion  coef¬ 
ficients  and  elastic  properties  do  not  exist.  Indeed,  for  multiphase  composites,  Schapery  [38] 
and  Rosen  and  Hashin  [35]  found  bounds  on  the  thermal  expansion  coefficients  in  terms 
of  the  stiffness  tensor.  Recently,  Gibiansky  and  Torquato  [11]  have  improved  upon  the 
Rosen-Hashin  bounds  using  the  so-called  translation  method.  This  improvement  was  actu¬ 
ally  motivated  by  the  topology  optimization  results  of  the  present  study.  To  our  knowledge, 
no  one  to  date  has  addressed  the  problem  of  systematically  designing  three-phase  materials 
with  extreme  isotropic  thermal  expansion  coefficients. 

In  this  paper,  we  show  how  composites  with  extremal  or  unusual  thermal  expansion  co¬ 
efficients  can  be  designed  using  a  three-phase  topology  optimization  method  based  on  the 
aforementioned  works  of  Sigmund.  The  three  phases  consist  of  two  different  material  phases 
and  void.  The  two  material  phases  can  have  different  elastic  and  thermal  expansion  coeffi¬ 
cients,  described  by  their  elastic  tensors  and  Cjjl  and  their  thermal  strain  coefficient 
tensors  o:-]^  and  ajp,  respectively.  The  basic  goal  is  to  maximize  or  minimize  components, 
or  combination  of  components,  of  the  effective  thermal  strain  tensor  «(•'  or  effective  stress 
tensor  subject  to  constraints  on  phase  volume  fractions,  effective  stiffness 

and  elastic  symmetry.  To  check  the  validity  of  the  optimization  procedure,  our  results  will 
be  compared  with  the  available  bounds  on  thermoelastic  properties  for  three-phase  materials 
[11,  35,  38]. 
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The  topology  optimization  procedure  proposed  here,  essentially  follows  the  steps  of  con¬ 
ventional  topology  optimization  procedures.  The  design  problem  is  initialized  by  defining 
a  design  domain  discretized  by  a  number  of  finite  elements.  The  optimization  procedure 
then  consists  in  solving  a  sequence  of  finite-element  problems  followed  by  changes  in  density 
and  material  type  of  each  of  the  finite  elements,  dependent  on  the  local  strain  energies.  For 
simple  compliance  optimization,  this  corresponds  to  adding  material  where  the  strain  energy 
density  is  high  and  removing  material  where  the  strain  energy  density  is  low. 

The  proposed  procedure  differs  from  the  conventional  approach  in  an  important  aspect. 
In  the  original  works  on  topology  optimization  (e.g.  Bendspe  and  Kikuchi  [4]  and  Kohn 
and  Strang  [23]),  elimination  of  ill-conditioning  and  the  existence  of  materials  in  elements 
of  intermediate  densities  was  ensured  by  using  so-called  ranked  materials  made  of  micro¬ 
scopically  oscillating  material  on  differing  length  scales  or  by  using  microstructures  with 
rectangular  holes.  Using  homogenization  methods  to  determine  the  effective  properties  of 
ranked  (or  microstructured)  materials  and  substituting  the  homogenized  parameters  into 
the  “macroscopic”  topology  optimization  problem,  lead  to  the  “homogenization  approach  to 
topology  optimization”.  In  this  paper,  however,  we  will  use  an  “artificial”  material  model 
for  intermediate  densities  which  means  that  the  material  properties  in  a  given  element  is 
simply  some  fraction  times  the  material  properties  of  solid  material.  As  long  as  we  end  up 
having  entirely  solid  material  or  void  in  each  element,  this  is  a  perfectly  valid  approach.  By 
using  the  “artificial-material”  model,  we  simplify  the  whole  design  procedure  significantly 
because  the  local  problems  of  determining  lamination  parameters  and  orientations  at  differ¬ 
ent  length  scales  of  the  ranked  materials  are  eliminated.  The  “artificial-material”  model  has 
been  used  by  several  authors  (e.g.  Rozvany  et  al.  [36],  Mlejnek  and  Schirrmacher  [29]  and 
Sigmund  [41]). 

At  each  step  of  the  topology  optimization  procedure,  we  have  to  determine  the  effective 
thermoelastic  properties  of  the  microstructure.  There  exist  several  methods  to  determine 
these  properties.  However,  because  the  topology  optimization  method  is  based  on  finite- 
element  discretizations,  and  because  the  finite-element  method  allows  easy  derivation  and 
evaluation  of  the  sensitivities  of  the  effective  properties  with  respect  to  design  changes,  we 
have  chosen  to  use  a  finite  element  based  numerical  homogenization  procedure  as  developed 
in  Bourgat  [6]  and  Guedes  and  Kikuchi  [14]. 

The  paper  is  organized  the  following  way.  In  section  2  we  describe  the  three-phase  topol¬ 
ogy  optimization  procedure  and  its  application  to  design  of  material  structures  with  extreme 
thermal  expansion.  The  sequential  linear  programming  method  used  to  solve  the  topology 
optimization  problem  is  described  in  section  2.2  and  numerical  implementation  issues  are 
discussed  in  section  2.3.  Solving  the  topology  optimization  procedure  as  proposed,  results  in 
solutions  with  finite-element  related  problems  such  as  checkerboards,  mesh-dependency  and 
local  optima.  These  problems  and  procedures  to  avoid  them  are  discussed  in  section  2.4. 
To  benchmark  the  optimization  procedure,  results  obtained  are  compared  to  bounds  for  two 
and  three  phase  materials;  the  available  bounds  are  listed  in  section  3.  The  performance 
of  the  design  procedure  is  demonstrated  by  several  examples  in  section  4.  The  calculation 
of  the  effective  thermoelastic  properties  using  a  numerical  homogenization  procedure  based 
on  the  finite-element  method  is  briefly  described  in  Appendix  A  and  the  sensitivity  analysis 
necessary  to  solve  the  design  problem  is  listed  in  Appendix  B. 
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2  Procedures  for  three-phase  topology  optimization 

This  section  describes  a  numerical  procedure  for  topology  optimization  of  three-phase  ma¬ 
terial  structures  in  two  dimensions.  A  sequential  linear  programming  problem  is  formulated 
to  solve  the  optimization  problem.  The  procedure  is  applied  to  the  design  of  material  struc¬ 
tures  with  extreme  thermal  expansion  properties.  For  the  derivations,  we  will  assume  that 
we  can  find  the  effective  stiffness  and  thermal  expansion  tensors  0-*^,  aj*^  and  using 
the  numerical  homogenization  method  described  in  Appendix  A.  The  sensitivity  analysis 
necessary  for  solving  the  problem  is  derived  in  Appendix  B.  At  the  end  of  the  section,  we 
discuss  implementation  issues  as  well  as  some  numerical  difficulties  and  how  to  avoid  them. 

Assuming  two-dimensional  linear  elasticity  (i.e.  small  strains),  perfect  bonding  between 
the  material  phases,  uniform  temperature  distribution  and  constant  material  properties,  the 
thermoelastic  behavior  of  materials  can  be  described  by  the  constitutive  relations  given  as 

~  ~  CijklOlkl^T  =  Cijkl^kl  —  (f) 

where  Cijki,  crij,  cw,  Okh  fhe  elasticity,  stress,  strain,  thermal  strain,  and  thermal  stress 

tensors,  respectively,  and  AT  is  the  temperature  change.  We  refer  to  aki  as  the  “thermal 
strain  tensor”  (the  resulting  strain  of  a  material  which  is  allowed  to  expand  freely)  and  to 
Pij  as  the  “thermal  stress  tensor”  (the  stress  in  a  material  which  is  not  allowed  to  expand). 
For  the  three-phase  composite  of  interest,  the  constitutive  equation  Eq.  (1)  is  valid  on  a 
local  scale  (with  superscripts  (0),  (1),  and  (2)  appended  to  the  thermoelastic  properties, 
i.e.  Cjjl],  and  and  the  macroscopic  scale  (with  superscript  (*)  appended  to  the 
properties).  In  the  latter  case,  the  stresses  and  strains  are  averages  over  local  stresses  and 
strains,  respectively,  i.e. 


where  overbar  denotes  the  volume  average.  The  effective  thermoelastic  properties,  Cjjl, 
and  of  the  three  phase  composite  are  computed  using  a  numerical  homogenization 
method  as  described  in  Appendix  A. 

The  goal  of  this  work  is  to  optimize  components  or  combinations  of  components  of  the 
effective  thermal  tensors  or  /Sj*^  by  distributing,  in  a  clever  way,  given  amounts  of  two 
material  phases  and  void  within  the  design  domain  representing  a  base  cell  of  a  periodic 
material.  In  other  words,  we  want  to  design  microstructural  topologies  that  give  us  some 
desirable  overall  thermoelastic  properties.  As  will  be  seen  later,  materials  with  extreme 
thermal  expansion  tend  to  have  low  overall  stiffness.  Thus,  for  practical  applications,  one 
must  bound  the  effective  stiffness  or  bulk  moduli  from  below.  It  should  also  be  possible  to 
specify  elastic  symmetries  such  as  orthotropy,  square  symmetry  or  isotropy  of  the  resulting 
materials. 
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An  optimization  problem  including  these  features  can  be  written  as 
Minimize  :  Some  function  of  or 

Variables  :  Distribution  of  two  material  phases  and  void  in  the  base  cell, 

Subject  to  :  Constraints  on  volume  fractions, 

:  Orthotropy,  square  symmetry  or  isotropy  constraints,  (3) 

:  Lower  bound  constraints  on  stiffness, 

:  Bounds  on  design  variables. 

2.1  Formulation  of  the  optimization  problem 

This  subsection  discusses  the  individual  parts  of  the  optimization  problem  defined  in  Eq.  (3). 
Objective  function 

The  objective  function  can  be  any  combination  of  the  thermal  coefficients 

or  An  example  will  be  the  case  where  we  want  to  minimize  the  isotropic  thermal 

expansion,  i.e.  the  sum  of  the  thermal  strain  coefficients  in  the  horizontal  and  the  vertical 
directions.  In  that  case,  the  objective  function  will  be  =  a^n  +  where  subscripts 

11  and  22  define  horizontal  and  vertical  directions,  respectively.  As  another  example  we 
might  consider  the  maximization  of  the  thermal  stress  coefficient  in  the  vertical  direction. 
In  this  case,  the  objective  function  will  be  =  -^22  where  the  minus  sign  is  used  to 

convert  the  maximization  problem  into  a  minimization  problem. 

Design  variables  and  mixture  assumption 

Phase  1  material  has  the  stiffness  tensor  Cjjh  and  the  thermal  strain  coefficient  tensor 
and  similarly  phase  2  material  has  the  material  tensors  C-Jh  and  a-^.  The  stiffness  tensor  of 
the  “void”  phase  is  taken  as  a  small  number  Xmin  times  C-Jl,  respectively,  where  x^in  - 
for  reasons  which  will  be  explained  later. 

The  material  type,  that  is,  material  phase  1,  phase  2  or  void,  can  vary  from  finite  element 
to  finite  element  as  seen  in  Fig.  1.  With  a  fine  finite-element  discretization,  this  allows  us  to 
define  complicated  bimaterial  topologies  within  the  design  domain.  Having  discretized  the 
design  domain  (the  periodic  base  cell)  with  N  finite  elements,  the  design  problem  consists 
in  assigning  either  phase  1,  2  or  void  to  each  element  such  that  the  objective  function  is 
minimized. 

Even  for  a  small  number  of  elements,  this  integer-type  optimization  problems  becomes 
a  huge  combinatorial  problem  which  is  impossible  to  solve.  For  a  small  design  problem 
with  N  =  100,  the  number  of  different  distributions  of  the  three  material  phases  would  be 
astronomical  (3^°*^  =  5  •  10^^).  As  each  function  evaluation  requires  a  full  finite  element 
analysis,  it  is  hopeless  to  solve  the  optimization  problem  using  random  search  methods 
such  as,  genetic  algorithms  or  simulated  annealing  methods,  which  use  a  large  number  of 
function  evaluations  and  do  not  make  use  of  sensitivity  information.  Following  the  idea  of 
standard  topology  optimization  procedures,  the  problem  is  therefore  relaxed  by  allowing  the 
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Design  domain  (base  cell) 


H  Phase  1  material 
^  Phase  2  material 

n  Void 


Periodic  material  structure 


Figure  1:  Design  domain  and  discretization  for  the  three-phase  topology  optimization  prob¬ 
lem.  Each  square  represents  one  finite  element  which  can  consist  of  either  phase  1  or  2 
material  or  void. 

material  at  a  given  point  to  be  a  mixture  of  the  three  phases.  This  makes  it  possible  to  find 
sensitivities  with  respect  to  design  changes,  which  in  turn  allows  us  to  use  mathematical 
programming  methods  to  solve  the  optimization  problem.  At  the  end  of  the  optimization 
procedure  however,  we  hope  to  have  a  design  where  each  element  is  either  void,  phase  1  or 
phase  2  material. 

Using  a  simple  artificicil  mixture  assumption,  the  local  stiffness  and  thermal  strain  coeffi¬ 
cient  tensors  in  element  e  can  be  written  as  a  function  of  the  two  design  variables  xl  and 


where  77  is  a  penalization  factor  discussed  later.  The  variable  x\  €  [x^in,  1]  can  be  seen  as  a 
local  density  variable  with  xl  =  Xmin  meaning  that  the  given  element  is  “void”  and  a;®  =  1 
meaning  that  the  given  element  is  solid  material.  The  variable  x^  €  [0, 1]  is  a  “mixture 
coefficient”  with  a;!  =  0  meaning  that  the  given  element  is  pure  phase  1  material  and  a:^  =  1 
meaning  that  it  is  pure  phase  2  material.  The  local  thermal  strain  tensor  alj{x2)  is  not 
dependent  on  the  density  variable  xj.  This  can  be  explained  by  the  fact  that  once  we  have 
chosen  the  local  material  mixture  (i.e.  the  value  of  x|),  the  thermal  strain  coefficient  does 
not  change  with  density. 

It  should  be  emphcisized  that  the  local  material  assumptions  Eqs.  (4)  only  are  valid  for 
the  design  variables  taking  the  extreme  values.^  Nevertheless,  during  the  design  process 
we  allow  intermediate  values  meaning  that  we  are  working  with  artificial  (non-existing) 

material  with  constitutive  behavior  close  to  the  described  could  be  realized  as  an  isotropic  porous 
triangular  microstructure 
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materials.  This  violation  is  not  critical  as  long  as  we  end  up  with  a  discrete  design  as 
discussed  in  the  introduction. 

Experience  shows  that  the  penalty  parameter  7}  should  be  given  values  ranging  from  3 
to  10  depending  on  the  design  problem.  The  influence  of  the  penalty  parameter  can  be 
explained  as  follows:  Let  us  assume  that  x\  =  ^  m  element  e.  The  local  stiffness  tensors 
dependence  of  [Eq.  (4)]  can  then  be  written  as  By  specifying  a 

value  of  rj  higher  than  one,  the  local  stiffness  for  fixed  <  1  is  lowered,  thus  making  it 
“uneconomical”  to  have  intermediate  densities  in  the  optimal  design. 


Constraints  on  volume  fractions 

Having  defined  the  design  variables  Xi  and  X2  above,  and  assuming  that  the  design  domain 
has  been  disrectized  by  N  finite  elements  of  volume  the  volume  fractions  of  the  three 
phases  can  be  calculated  as  the  sums 


^  e=\ 


eO)  =  1  _  c(')  - 


where  Y  is  the  volume  of  the  base  cell.  For  a  specific  design  problem,  we  might  want  to 
constrain  the  volume  fractions  of  the  phases.  This  can  be  done  by  defining  two  volume 
fraction  constraints  as 


,  ^  ^maxi  ^min  ~  ^  —  ''mon 


where  c^L,  c^nf  4}llx  ^mlx  are  lower  and  upper  bounds  on  the  volume  fractions  of 
material  1  and  2,  respectively.  By  setting  the  lower  bound  constraint  equal  to  the  upper 
bound  constraint,  it  is  possible  to  fix  the  volume  fractions  of  the  individual  phases. 


Isotropy  or  square  symmetry  constraints 

For  the  purpose  of  designing  materials  with  either  orthotropic,  square  symmetric  or  isotropic 
elastic  parameters,  such  constraints  must  be  implemented  in  the  optimization  problem.  Or- 
thotropy  of  the  materials  can  be  obtained  simply  by  specifying  at  least  one  geometrical 
symmetry  axis  in  the  base  cell.  Assuming  that  a  material  structure  is  orthotropic,  the  condi¬ 
tion  for  square  symmetry  of  the  elasticity  tensor  is  that  Cim  —  C2222  —  0?  conditions 

for  isotropy  of  the  elasticity  tensor  under  plane  stress  assumption  are  that  Cmi  —  C2222  =  0 
and  (Cllli  +  CS22)  -  +  2C'J;}2)  -  0.  Finally,  the  condition  for  thermal  expansion 

isotropy  is  that  =  0  and  0:^2^  =  0.  These  conditions  are  difficult  to  implement 

as  equality  constraints  in  an  optimization  problem  because  the  starting  guess  might  be  in¬ 
feasible  (i.e.  anisotropic).  Therefore,  it  is  chosen  to  implement  the  constraints  as  a  penalty 
function  added  to  the  cost  function.  The  penalty  function  is  defined  as  the  squared  error  in 
obtaining  either  square  symmetry,  elastic  or  thermal  isotropy,  times  the  penalization  factors 
n,  rz  and  rg,  respectively.  It  should  be  noted  here  that  three  60  degree  symmetry  lines  of  a 
microstructure  is  a  sufficient  but  not  a  necessary  condition  for  isotropy.  Indeed,  this  paper 
shows  examples  of  isotropic  material  structures  with  only  one  line  of  symmetry. 


The  errors  in  obtaining  square  symmetry  or  isotropy,  respectively,  can  be  written  as 
Error  -  ~  *^2222)^ 

p  [(Cjni  +  C2222)  ~  2(f?n22  +  2C1212)] 

E’-ror<„  =  -  - <-  + Error.,. 

I'-'llll  +  ^2222! 


Expressions  similar  to  Error and  Error^so  are  also  known  in  the  literature  of  composite 
materials  (e.g.  Christensen  [7])  as  the  practical  composite  parameters  U2  and  f/a. 

The  error  in  obtaining  isotropic  thermal  expansion  can  be  defined  as 


Errori}2^j.jji  — 


)^  +  (4V-4 


Lower  bound  constraints  on  effective  stiffness 

As  will  be  seen  later,  extreme  thermal  properties  can  be  obtained  if  we  allow  the  overall 
stiffness  of  the  material  to  be  small.  Low  stiffness  is  generally  undesirable  and  therefore  we 
will  introduce  a  lower  bound  constraint  on  the  the  directional  Young’s  moduli  E[*^  and  E^*^ 
or  on  the  bulk  modulus  of  the  material. 

Such  constraints  can  be  written  as  gi^-min)  <  9i{Cjjli).  For  isotropic  materials,  we  have 
a  lower  bound  constraint  on  the  bulk  modulus  +  CjiLl- 

For  anisotropic  materials^  we  might  want  to  constrain  the  value  of  the  horizontal  or  vertical 
Young  s  moduli,  i.e.  <  -£^1  ^  —  (C'ii22)^/C'22L]  or  <  ^^2*  =  ~ 

Lower  bound  constraints  on  design  variables 

For  computational  reasons  (singularity  of  the  stiffness  matrix  in  the  finite  element  formula¬ 
tion),  the  lower  bound  on  design  variable  is  set  to  x^in]  not  zero  {xmin  =  10"'*).  Numerical 
experiments  show  that  the  “void”  regions  have  practically  no  structural  significance  and  can 
be  regarded  as  real  void  regions.  The  bounds  on  the  design  variables  can  thus  be  written  as 
0  <  Xmin  <  <  1  and  0  <  a:|  <  1. 

The  final  optimization  problem 

An  optimization  problem  including  above  mentioned  features  can  now  be  written  as 
Minimize  :  $(xi,  X2)  =  /(q:S*\  /?|;^)  +  r^Error.q,  +  r^Erran^o  +  rzErrortherm, 
subject  to  :  p, <  ^,((7,^*^),  i  =  1, . . . ,  M, 

•  ^mtn  —  ^  —  ^maxi  (,"1 

•  <  c(2)  <  c(2) 

•  ^min  —  ^  —  ^maxi 

•  0  ^min  ^  5Ci  ^  1, 

:  0<X2<1, 
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where  Xi  and  X2  are  the  iV-vectors  containing  the  design  variables  and  the  three  penalty 
parameters  r,-  can  be  set  to  zero  or  non-zero  values,  depending  on  the  desired  isotropy  type. 


2.2  Sequential  linear  programming  method 

Topology  optimization  problems  in  the  literature  often  consist  in  the  optimization  of  a  simple 
energy  functional  (e.g.  compliance  or  eigenfrequencies)  with  a  single  constraint  on  material 
resource,  and  these  problems  can  therefore  be  solved  very  efficiently  using  optimality  criteria 
methods.  In  this  paper,  however,  we  are  considering  several  different  objective  functions  and 
multiple  constraints  which  can  not  be  written  in  energy  forms  and  therefore  it  will  be  cum¬ 
bersome  if  not  impossible  to  formulate  the  optimization  problem  as  an  optimality  criteria 
problem.  Instead  we  will  use  a  mathematical  programming  method  called  sequential  linear 
programming  (SLP),  which  consists  in  the  sequential  solving  of  an  approximate  linear  sub¬ 
problem,  obtained  by  writing  linear  Taylor  series  expansions  for  the  objective  and  constraint 
functions.  The  SLP  method  was  successfully  used  in  optimization  of  truss  structures  by 
Pedersen  [32]  and  was  evaluated  as  a  robust,  efficient  and  easy  to  use  optimization  algorithm 
in  a  review  paper  by  Schittkowski  [39]. 

Using  the  sequential  linear  programming  method,  the  optimization  problem  Eq.  (9)  is 
solved  iteratively.  In  eaoh  iteration  step,  the  optimization  problem  is  linearized  around  the 
current  design  point  {xi,X2}  using  the  first  part  of  a  Taylor  series  expansion  and  the  vector 
of  optimal  design  changes  {Axi,  AX2}  is  found  by  solving  the  linear  programming  problem 


f  5$  r 

Minimize  :  {^XijAxj}, 

subject  to  :  gi^min)  -  9i  <  AX2},  i  =  1, . . . ,  M, 


:  {Axil,Ax2l}  <  {Axi,Ax2}  <  {Axjc/,  Ax2y}, 


where  Axil,  Ax2Z„  Axit;  and  Ax2r7  are  move-limits  on  the  design  variables.  The  move-limits 
are  adjusted  for  the  absolute  limits  given  in  Eqs.  (9). 

The  applied  move-limit  strategy  is  important  for  the  stable  convergence  of  the  algorithm. 
Here  we  use  the  simple  rule  that  the  move-limit  for  a  specific  design  variable  is  increased 
by  a  factor  of  1.4  if  the  change  in  the  design  variable  has  the  same  sign  for  two  subsequent 
steps.  Similarly  the  move-limit  is  decreased  by  a  factor  of  0.6  if  the  change  in  the  design 
variable  has  opposite  signs  for  two  subsequent  steps. 

The  sensitivities,  which  are  necessary  to  solve  the  linearized  sub-problem  Eq.  (10),  are 
derived  in  Appendix  B,  and  are  calculated  locally  for  each  element.  This  means  that  no 
additional  finite-element  problems  have  to  be  solved  to  find  the  sensitivities  needed. 
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Figure  2:  Flowchart  of  the  design  algorithm 

2.3  Numerical  implementation  issues 

This  subsection  describes  the  numerical  implementation  of  the  three-phase  topology  op¬ 
timization  problem  including  the  finite-element  discretization  and  procedures,  the  linear 
programming  package  DSPLP  [16]  from  the  SLATEC  library,  control  of  move-limits  and  a 
flowchart  of  the  procedure. 

A  flowchart  of  the  design  algorithm  is  shown  in  Fig.  2.  The  individual  steps  of  the  design 
procedure  are  described  in  the  following. 

Initialization 

First,  we  initialize  the  design  problem  by  selecting  the  objective  function,  specifying  a  lower 
bound  on  the  stiffness,  selecting  isotropy  type  and  symmetry  lines.  We  also  choose  the  design 
domain  discretization,  using  900  or  3600  4-node  linear  displacement  finite  elements,  corre¬ 
sponding  to  30  by  30  or  60  by  60  element  discretizations,  depending  on  accuracy  demands 
and  available  computing  time.  To  save  computer  time,  a  design  problem  can  first  be  solved 
on  a  30  by  30  element  mesh.  When  a  solution  has  been  reached,  each  of  the  elements  are 
divided  into  four  and  the  procedure  is  continued  until  convergence. 

Starting  guess 

Starting  distribution  of  densities  and  material  types  (i.e.  starting  values  of  the  design  variable 
vectors  Xi  and  X2)  is  up  to  the  user.  Having  absolutely  no  idea  of  what  the  solution  will  look 
like,  a  random  distribution  of  densities  and  material  types  is  chosen  as  the  starting  guess.  If 
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the  user  has  an  idea  of  what  the  solution  will  look  like  or  he  has  an  old  solution  to  a  similar 
problem,  considerable  amount  of  computing  time  is  saved  by  using  this  (old)  topology  as  a 
starting  guess. 

Homogenization  step 

The  equilibrium  equations  Eqs.  (17)  for  the  homogenization  problem  derived  in  Appendix  A, 
are  solved  using  the  finite-element  method  applied  to  calculation  of  effective  material  prop¬ 
erties  in  Bourgat  [6]  and  Guedes  and  Kikuchi  [14]. 

Sensitivity  analysis 

The  sensitivity  analysis  necessary  to  solve  the  linear  programming  problem  in  Eq.  (10)  is 
derived  in  Appendix  B.  The  computation  of  the  sensitivities  is  fast  because  they  can  be 
found  from  the  strain  fields  already  computed  by  the  homogenization  procedure. 

Linear  programming  problem 

The  linear  programming  problem  Eq.  (10)  is  solved  using  a  linear  programming  solver 
DSPLP  [16]  from  the  SLATEC  library.  As  the  optimization  is  non-sparse,  the  DSPLP 
routine  is  invoked  with  an  option  for  no  sparsity.  Nevertheless,  the  routine  has  proven  faster 
and  demands  less  storage  space  than  other  LP-algorithms  tests. 

Convergence 

The  iterative  design  procedure  is  repeated  until  the  change  in  each  design  variable  from  step 
to  step  is  lower  than  10"'*  (by  experience). 

2.4  Problems  related  to  topology  optimization 

This  subsection  discusses  some  numerical  difficulties  due  to  the  finite- element  discretization, 
namely  checkerboard  patterns,  mesh-dependencies  and  local  minima. 

The  checkerboard  and  mesh-dependency  problems 

Applying  the  topology  optimization  method  to  different  design  problems,  one  often  encoun¬ 
ters  regions  of  alternating  solid  and  void  elements,  referred  to  as  checkerboards,  in  the 
“optimal  solutions”.  The  regions  are  seen  in  many  works  on  general  topology  optimization 
and  it  was  earlier  believed,  that  such  regions  represented  optimal  microstructure  on  the 
finite-element  level.  However,  two  recent  papers  by  Jog  and  Haber  [21]  and  Diaz  and  Sig¬ 
mund  [9]  conclude,  that  regions  with  checkerboard  patterns  have  artificially  high  (numerical) 
stiffness  (higher  than  the  theoretical  bounds)  and  can  be  explained  by  poor  numerical  mod¬ 
elling  of  the  stiffness  of  checkerboards  by  lower  order  finite  elements.  Both  papers  conclude, 
that  checkerboards  are  certainly  prone  to  appear  in  topology  optimization  using  four-node 
finite  elements,  as  here,  but  also  using  higher  order  elements  such  as  nine-node  quadratic 
displacement  finite  elements. 
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Another  problem,  due  to  the  finite-element  discretization,  is  mesh-dependency,  which 
refers  to  the  non-convergence  of  solutions  with  mesh-refinement.  Refining  the  finite- element 
mesh  should  ideally  result  in  the  same  topology  as  for  a  coarse  mesh  but  with  better  defini¬ 
tions  of  the  boundaries  between  the  material  phases.  However,  a  refinement  does  result  in  a 
solution  with  a  more  complicated  (finer)  microstructure.  Methods  to  avoid  this  problem  have 
been  suggested  in  three  recent  papers.  Jog,  Haber  and  Bends0e  [15]  suggest  to  introduce  a 
constraint  on  global  perimeter.  In  a  paper  on  bone  remodelling,  closely  related  to  topology 
optimization  methods,  Mullender,  Huiskes  and  Weihnans  [30]  suggest  a  mesh-independency 
algorithm  that  assumes  that  bone  growth  at  a  point  is  dependent  on  the  loads  in  a  (mesh- 
independent)  neighborhood  of  the  point.  A  method  related  to  the  latter  approach,  but  with 
different  origin,  is  proposed  in  Sigmund  [41],  who  performs  the  density  update  based  on  low- 
pass  filtered  strain  energy  fields.  To  avoid  the  checkerboard  and  mesh- dependency  problems, 
we  use  the  method  suggested  in  Sigmund  [41] 

Local  minima 

The  topology  optimization  problem  is  very  prone  to  converge  to  local  minima.  However, 
introducing  the  mesh-independency  algorithm  (Sigmund  [41])  makes  it  possible  to  prevent 
this  problem  to  a  certain  extent.  Solving  a  cell  design  problem  is  typically  done  as  follows. 
First  we  solve  the  optimization  problem  with  a  low  value  of  the  low-pass  filter  parameter, 
i.e.,  we  do  not  allow  rapid  variation  in  the  element  densities.  This  results  in  a  design  with 
large  areas  of  intermediate  densities  but  it  also  prevents  the  design  in  converging  to  a  local 
minimum  (binary  design).  Gradually,  we  increase  the  low-pass  filter  parameter,  in  turn 
letting  the  design  problem  converge.  In  that  way,  we  gradually  arrive  at  a  solution  which  is 
entirely  binary  and  which  is,  hopefully,  a  global  optimum.  To  make  sure  that  the  actually 
obtained  microstructures  are  global  optima  indeed,  the  same  optimization  problem  is  always 
solved  using  differing  starting  guesses,  move-limit  strategies  and  choices  of  low-pass  filter 
parameter  and  penalty  parameter  t).  However,  as  will  be  seen  later,  topologically  different 
solutions  with  similar  values  of  the  objective  function  have  been  found  when  solving  specific 
design  problems.  Solutions  which  are  ’’shifted”  (translated  half  a  base  cell  dimension)  of 
other  solutions  have  also  been  encountered.  The  fact  that  the  effective  properties  of  the 
design  examples  are  close  to  theoretical  bounds  supports  our  belief,  that  we  are  finding  the 
optimal  topologies  with  the  proposed  design  procedure. 

Computing  time 

One  design  iteration,  typically  takes  three  seconds  (30  by  30  element  discretization)  and  20 
seconds  (60  by  60  element  discretization)  on  an  Indigo  2  work  station.  To  arrive  at  an  opti¬ 
mal  solution,  depending  on  starting  guess,  several  thousand  iterations  are  needed.  Including 
interaction  by  the  user,  a  full  design  process  may  take  two  working  days.  Comparing  the 
computing  time  and  number  of  iterations  with  other  works  in  topology  optimization  it  can 
be  concluded  that  the  present  procedure  is  extremely  slow.  However,  these  other  works  in 
topology  optimization  usually  consider  self-adjoint  loading  problems  and  statically  determi¬ 
nant  structures  where  a  solution  can  be  found  often  within  a  few  design  iterations.  This  is 
due  to  the  fact  that  the  stress  fields  of  statically  determinate  structures  do  not  change  much 
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with  design.  In  this  paper,  however,  the  optimal  strain  and  stress  fields  for  the  given  design 
problem  are  unknown  in  the  beginning  and  must  emerge  gradually  during  the  design  process 
together  with  the  optimal  material  distribution. 


3  Rigorous  expressions  and  bounds  on  effective  ther¬ 
moelastic  properties 

Rigorous  expressions  for  the  effective  thermal  expansion  coefficients  of  two-phase,  isotropic 
composites  and  rigorous  bounds  on  the  effective  coefficients  of  three-phase,  isotropic  com¬ 
posites  will  serve  to  benchmark  the  design  algorithm.  For  simplicity,  we  assume  that  the 
constituent  phases  are  isotropic  which  implies  that  they  can  be  described  by  their  Young’s 
moduli  and  E^'^\  their  Poisson’s  ratios  and  and  their  thermal  strain 

coefficients  and  It  is  also  assumed  that  the  composite  is  macroscopicically 

isotropic.  The  bulk  and  shear  moduli  of  the  phases  are  then 

i  =  0,1,2.  (11) 


2(1 -I/O))’ 


2(1  -l-t'W)’ 


Two-phase  materials 

The  effective  thermal  strain  coefficient  of  a  two-phase  isotropic  material  is  explicitly 
given  in  terms  of  the  effective  bulk  modulus  (Levin  [27]  or  Rosen  and  Hashin  [35]) 

-  A:(*))  -  al2)jfc(2)  (fclD  _ 

"  fc(*)  (fc(2)  -  fc(i))  ■ 

The  best  bounds  on  the  isotropic  effective  bulk  modulus  k^*\  given  volume  fraction  infor¬ 
mation  only,  were  derived  by  Hashin  and  Shtrikman  [17]  and  read 


where  bars  denote  averaged  values  and  pmin  and  pmax  are  the  minimum  and  maximum  shear 
moduli  [Eq.  (11)]  of  the  three  phases,  respectively.  Bounds  for  the  thermal  strain  coefficient 
are  therefore  obtained  by  inserting  the  upper  and  lower  bounds  for  the  bulk  modulus  Eq.  (13) 
into  Eq.  (12). 

Three-phase  materials 

Bounds  on  the  effective  thermal  strain  coefficient  of  three-phase,  isotropic 

composites  were  found  by  Schapery  [38]  and  Rosen  and  Hashin  [35]  and  read 
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(14) 


where 


Note  that  there  is  a  typographical  error  in  the  equation  [Eq.  (14)]  in  Rosen  and  Hashin  [35] 
as  well  as  in  the  corresponding  formula  in  the  text  by  Christensen  [7]. 

A  bounded  domain  of  possible  effective  bulk  moduli  and  thermal  strain  coefficients  for  a 
specific  choice  of  constituent  phases  is  shown  in  Fig.  4,  We  found  that  the  proposed  design 
method  did  not  yield  pairs  that  were  close  to  the  Schapery-Rosen-Hashin  bounds 

[Eq.  (14)].  There  were  two  possible  explanations  for  this  discrepancy:  either  the  design 
method  could  not  find  the  optimal  solutions,  or  the  bounds  themselves  could  be  improved 
upon.  Indeed,  the  latter  explanation  turned  out  to  be  true. 

Inspired  by  above  mentioned  discrepancy  Gibiansky  and  Torquato  [11]  recently  found 
improved  bounds,  which  are  also  shown  in  Fig.  4.  The  new  bounds  can  be  written  as 

(^(*)  -  ^^^^)(^^^^  +  j  ±  (15) 

where 


and  and  are  the  lower  and  upper  Hashin-Shtrikman  bounds  on  bulk  modulus  as 
given  in  Eq.  (13).  As  will  be  seen  in  the  subsequent  section,  the  solutions  obtained  by  the 
design  procedure  axe  very  close  to  the  new  bounds. 

Examination  of  the  thermoelastic  bounds  in  Fig.  4  reveals  that  extreme  values  (e.g. 
negative  values)  of  thermal  strain  coefficients  only  are  possible  for  low  bulk  moduli.  If  we 
simply  tried  to  rainimize/maximize  the  thermal  strain  we  would  end  up  with  a  very 
weak  material.  Therefore,  there  is  a  tradeoff  between  extremizing  thermal  strain  coefficients 
on  the  one  hand  and  ending  up  with  a  stiff  material  on  the  other.  This  problem  will  be 
discussed  in  more  detail  in  the  next  section. 
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4  Design  examples 

In  this  section,  we  will  first  discuss  design  examples  with  mixtures  of  hypothetical  materials. 
These  examples  are  used  to  benchmark  the  design  algorithm  for  two-  and  three-phase  design. 
We  will  also  study  other  design  examples  that  utilize  real  materials  as  constituent  phases. 

Plotting  results 

During  the  iterative  procedure,  a  Postscript  plot  of  the  topology  is  generated  every  10  iter¬ 
ations.  The  plot  shows  the  current  density  and  material  distribution  in  the  base  cell,  thus 
allowing  the  user  to  follow  the  evolution  of  the  microstructure  and  interact  if  necessary.  The 
plots  in  the  following  sections  show  the  optimal  density  and  material  type  distributions  for 
the  different  design  problems.  If  an  element  is  dominating  material  phase  two  (i.e.,  >  0.5), 

the  element  is  illustrated  by  a  cross  with  grey  scale  denoting  the  density  white  means 
void  (^5  =  x,„in)  and  black  means  solid  (xf  =  1).  If  the  element  is  dominating  material 
phase  one  (i.e.,  <  0.5),  it  is  shown  as  a  filled  rectangle  with  grey  values  interpreted  as 

before.  For  all  examples,  we  both  show  the  resulting  topologies  represented  by  a  single  base 
cell  (the  design  domain)  and  as  a  repeated  microstructure  consisting  of  3  by  3  base  cells. 

4.1  Comparison  with  two-phase  bounds 

The  two-phase  bounds  given  by  Eqs.  (12)  and  (13)  are  used  for  benchmarking  the  design 
algorithm  for  two-phase  design.  The  non-dimensionalized  material  data  for  the  two  phases 
are  chosen  as  =  10,  =  0.3,  =  iQ  and  =  0.5. 

The  bounds  on  the  thermal  strain  coefficient  for  this  mixture  are  found  from  Eq.  (12) 
to  be  given  by  6.3412  <  <  8.2524  and  the  bounds  on  the  bulk  modulus  are  found 

from  Eq.  (13)  to  be  given  by  0.0996  <  <  0.1690. 

First  we  try  to  maximize  the  effective  thermal  strain  coefficient  which  through  Eq.  (12) 
corresponds  to  maximizing  the  effective  bulk  modulus.  Specifying  macroscopic  isotropy  and 
horizontal  and  vertical  geometric  symmetry,  the  attained  values  with  a  60  by  60  element 
discretization  are  =  8-2349  and  kl;lJk^^^  =  0.1680.  If  we  try  to  minimize  the 

thermal  strain  coefficient  we  just  get  the  “inverted”  microstructure,  meaning  that  the  two 
domains  are  interchanged.  The  actual  numbers  obtained  for  this  case  are  =  6.3907 

and  =  0.1007.  The  optimal  topology  of  the  microstructure  for  maximum  thermal 

strain  value  and  isotropy  constraint  is  shown  in  Fig.  3  (top). 

Relaxing  the  isotropy  requirement  by  only  specifying  square  symmetry  (but  still  spec¬ 
ifying  horizontal  and  vertical  geometric  symmetry),  the  attained  values  with  a  60  by  60 
element  discretization  are  =  8.2320  and  k^l^/k^^'>  =  0.1677,  =  6.3733  and 

=  0.1004,  respectively.  The  actual  topology  of  the  optimal  two-phase  microstruc¬ 
ture  is  shown  in  Fig.  3  (bottom).  As  expected,  the  optimal  microstructural  topologies 
resemble  the  energy  minimizing  microstructures  of  Vigdergauz  [45] . 

It  should  be  noted  that  we  get  the  same  (except  for  numerical  errors)  effective  properties 
for  the  isotropic  material  and  square  symmetric  microstructures.  This  means,  that  the 
optimal  rigidity  Vigdergauz  microstructures  can  be  made  isotropic  by  using  the  new  geometry 
in  Fig.  3  (top)  instead  of  the  one  in  Fig.  3  (bottom).  Not  surprisingly,  it  also  shows  that 
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Figure  3;  Optimal  microstructures  for  two-phase  design  problem.  Maximization  of  thermal 
strain  coefficient  with  and  without  macroscopic  isotropy  constraint  (top  and  bottom  respec¬ 
tively).  The  filled  regions  consist  of  low  expansion  material  (phase  1)  and  the  cross-hatched 
regions  consist  of  high  expansion  material  (phase  2). 

solutions  with  widely  different  topologies  can  have  the  same  values  of  the  objective  function. 


4.2  Comparison  with  three-phase  bounds 

The  three-phase  bounds  given  by  Eqs.  (13),  (14)  and  (15)  are  used  for  benchmarking  the 
design  algorithm  for  three-phase  design.  The  material  data  for  the  two  phases  are  chosen  as 
£'{i)/£'(2)  _  2^  j,(i)  _  j,(2)  _  Q  volume  fractions  are  prescribed  to 

be  =  0.25  (i.e.  =  0.5).  ^'ote  that  the  volume  fractions  Ci  are  held  fixed  for  this 

hypothetical  composite,  to  allow  for  comparison  with  the  bounds  and  for  easy  interpretation 
of  the  results. 


Effective  thermal  strain  coefficient  a 


Figure  4:  Bounds  for  three-phase  design  example.  The  circles  with  letters  a-d  denote  the 
obtained  values  for  the  microstructures  shown  in  Figs.  5  and  6. 


Four  three-phase  design  examples,  constrained  to  be  elastically  isotropic,  are  considered 
as  follows: 

(a)  Minimization  of  the  isotropic  thermal  strain  coefficient  with  a  lower  bound 

constraint  on  the  effective  bulk  modulus  given  as  10%  of  the  theoretically  attainable 

bulk  modulus,  i.e.  —  0.0258.  Horizontal  geometric  symmetry  is  specified. 

(b)  Same  as  design  example  (a)  but  with  horizontal,  vertical  and  diagonal  (geometric) 

symmetry. 

(c)  Maximization  of  bulk  modulus  k^*yk^^'>  for  fixed  zero  thermal  expansion  =  0. 

Horizontal  geometric  symmetry  is  specified. 

(d)  Maximization  of  isotropic  thermal  stress  coefficient  with  horizontal,  vertical 

and  diagonal  geometric  symmetry. 

The  old  and  new  theoretical  bounds  are  given  by  Eqs.  (13)  and  (15),  respectively,  and  they 
are  shown  Fig.  4.  In  examples  (a)  and  (b),  the  lower  bound  on  the  possible  thermal  strain 
coefficient  is  —5.567  <  In  example  (c),  the  upper  bound  on  possible  bulk  modulus 

for  zero  thermal  expansion  is  <  0.0692.  The  upper  bound  on  the  thermal  stress 

coefficient  in  design  example  (d)  is  <  3.15  (for  =  0.237). 

The  resulting  topologies  are  shown  in  Figs.  5  and  6  and  their  effective  properties  are 
shown  in  Table  1  and  plotted  as  small  circles  in  Fig.  4.  Studying  the  graph  in  Fig.  4,  we 
see  that  the  obtained  effective  values  are  far  away  from  the  original  Schapery-Rosen-Hashin 
bounds.  This  discrepancy  inspired  Gibiansky  and  Torquato  to  try  to  improve  the  bounds 
and  indeed  improvement  was  possible  as  seen  in  Fig.  4.  The  effective  values  of  the  examples 
(a)-(d)  are  still  somewhat  away  from  the  improved  bounds.  This  can  be  explained  by  the 
fact  that  the  new  bounds  by  Gibiansky  and  Torquato  have  not  been  proven  to  be  optimal. 
Furthermore,  it  is  our  experience  that  a  finer  finite-element  mesh  makes  it  possible  to  get 
closer  to  the  bounds.  In  example  (a),  the  minimum  thermal  strain  coefficient  obtained  for 
a  30  by  30  mesh  is  =  —3.59  and  =  —4.17  for  the  60  by  60  element  discretization 
shown  in  Fig.  5.  Due  to  computer  time  limitations,  it  has  not  been  possible  to  try  out  finer 
discretizations. 

The  actual  mechanisms  behind  the  extreme  thermal  expansion  coefficients  of  the  ma- 
terialstructures  can  be  difficult  to  understand.  To  visualize  one  of  the  mechanisms,  the 
(exaggerated)  displacements,  due  to  an  increase  in  temperature  of  the  microstructure  in 
Fig.  5  (bottom),  is  shown  in  Fig.  7.  Studying  Fig.  7,  we  note  that  there  appears  to  be 
contact  between  parts  of  the  microstructure.  This  contact  is  only  due  to  the  magnification 
of  the  displacements  used  in  the  illustration.  The  simple  linear  modelling  used  here  can  not 
take  such  problems  into  account.  Nevertheless,  it  would  be  interesting  to  extend  the  analysis 
to  include  non-linear  behavior  including  contact,  which  would  open  up  for  a  whole  new  world 
of  interesting  design  possibilities.  We  will  leave  these  extensions  to  future  studies. 

When  allowing  low  bulk  moduli  [as  in  examples  (a)  and  (b)],  the  main  mechanism  behind 
the  extreme  (negative)  thermal  expansion  is  the  reentrant  cell  structure  having  bimaterial 
components  which  bend  and  cause  large  deformation  when  heated.  The  bimaterial  interfaces 
of  design  examples  (a)  and  (b)  bend  and  make  the  cell  contract,  similar  to  the  behavior 
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Figure  5:  Examples  (a)  (top)  and  (b)  (bottom):  Optimal  microstructures  for  minimization 
of  effective  thermal  strain  coefficient  corresponding  to  the  circles  a  and  6  in  Fig.  4,  respec¬ 
tively.  The  white  regions  denote  void  (phase  0),  the  filled  regions  consist  of  low  expansion 
material  (phase  1)  and  the  cross-hatched  regions  consist  of  high  expansion  material  (phase 


Figure  6:  Examples  (c)  (top)  and  (d)  (bottom):  Optimal  microstructures  for  maximization 
of  bulk  modulus  with  zero  thermal  expansion  (top)  and  maximization  of  effective  thermal 
stress  coefficient  (bottom)  corresponding  to  the  circles  c  and  d  in  Fig.  4,  respectively.  The 
white  regions  denote  void  (phase  0),  the  filled  regions  consist  of  low  expansion  material 
(phase  1)  and  the  cross-hatched  regions  consist  of  high  expansion  material  (phase  2). 


Figure  7:  Thermal  displacement  of  microstructure  in  Fig.  5  (bottom). 


Example 

Objective/ 

Constraint 

(bound) 

mu 

BUSH 

(a),  Fig.  5 

Min.  a''*^/a''^'^ 

>  0.0258 

0.0258 

0.039 

-4.17 

(-5.567) 

(b),  Fig.  5 

Min. 

>  0.0258 

0.0258 

0.51 

-4.02 

(-5.567) 

(c),  Fig.  6 

Max. 

ja^^^  <  0.0 

0.0692 

(0.0814) 

0.54 

0 

(d),  Fig.  6 

Max.  ^ 

0.243 

0.51 

3.01 

(3.15) 

Table  1:  Thermoelastic  parameters  for  optimal  three-phase  microstructures  composed  of 
hypothetical  materials  compared  with  the  bounds.  The  white  regions  denote  void,  the  filled 
regions  consist  of  low  expansion  material  (phase  1)  and  the  cross-hatched  regions  consist  of 
high  expansion  material  (phase  2). 


of  negative  Poisson’s  ratio  materials  [24].  If  a  higher  effective  bulk  modulus  is  specified,  as 
in  example  (c),  the  intricate  bi-material  mechanisms  are  less  pronounced  resulting  in  a  less 
extreme  expansion  (a,  =  0).  Finally,  maximizing  the  expansive  stress,  as  in  example  (d), 
results  in  a  structure  without  bimaterial  mechanisms,  where  the  high  expansion  phase  (cross 
hatched  phase)  is  arranged  such  that  it  maximizes  the  horizontal  and  vertical  expansion. 

Design  examples  (a)  and  (b)  in  Fig.  5  demonstrate  how  two,  topologically,  very  differ¬ 
ent  microstructures  can  have  (almost)  the  same  value  of  the  objective  function.  The  only 
difference  between  the  two  examples  is  the  specified  geometric  symmetry. 

In  order  to  check  the  validity  of  the  effective  values  computed  by  the  homogenization  code, 
we  will  make  a  comparison  between  a  simple  example  and  an  idealized  beam  model.  Figure  8 
shows  the  optimal  microstructure  (e)  for  maximization  of  the  thermal  strain  coefficient  in 
the  vertical  direction  022^  We  specify  geometrical  symmetry  about  horizontal  and  vertical 
axes,  effective  vertical  stiffness  component  C2222  ^  0-07  and  volume  fractions  and 
less  or  equal  0.2.  The  obtained  values  are  =  —20.7  and  0:22^  =  33.5.  The  actual  material 
volume  fractions  in  the  optimal  topology  are  c^^l  =  0.167  and  =  0.200  and  effective 
horizontal  and  vertical  stiffness  components  are  =  0.070,  respectively. 

To  check  the  obtained  values,  we  calculate  the  effective  thermal  strain  tensor  using  an 
idealized  beam  model  in  Appendix  C.  The  predicted  thermal  strain  coefficients  and  bulk 
modulus,  are  =  -21,  ay  =  32  and  =  0.085,  respectively.  Good  agreement  between 
the  idealized  beam  model  and  the  effective  properties  is  observed. 
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Figure  8:  Example  (e):  Optimal  microstructure  for  maximization  of  thermal  strain  in  the 
vertical  direction  oc^22-  The  white  regions  denote  void  (phase  0),  the  filled  regions  consist 
of  low  expansion  material  (phase  1)  and  the  cross-hatched  regions  consist  of  high  expansion 
material  (phase  2). 

4.3  Mixtures  of  real  materials 

For  the  design  of  new  materials  with  extreme  thermal  expansion  coefficients,  the  two  base 
materials  should  be  of  equal  stiffness  but  widely  differing  thermal  strain  coefficients.  Two 
materials  fulfilling  this  requirement  are  isotropic  Invar  (Fe-36%Ni)  and  Nickel  as  discussed 
in  the  introduction.  For  the  next  design  examples,  the  volume  fractions  of  the  material  phases 
are  unconstrained.  This  will  allow  for  a  wider  range  of  minimum  and  maximum  values,  in 
contrast  to  the  hypothetical  examples  (a)-(e)  in  which  the  volume  fractions  were  fixed. 

The  material  properties  of  Invar  and  Nickel  can  be  found  in  the  ASM-Handbook  [1].  The 
Young’s  moduli  are  150  GPa  and  200  GPa,  respectively,  Poisson’s  ratios  are  0.31  for  both, 
and  the  thermal  expansion  coefficients  are  0.8^m/(mA")  and  13.4/im/(mAr),  respectively. 

(f)  Minimization  of  the  isotropic  thermal  stress  coefficient  Horizontal  geometric  sym¬ 

metry  is  specified. 

(g)  Minimization  of  the  vertical  thermal  stress  coefficient  ^22  ■  Horizontal  and  vertical  sym¬ 

metry  is  specified. 

(h)  Minimization  of  the  vertical  thermal  stress  E2  •>4?.  Horizontal  and  vertical  symmetry 

is  specified. 

(i)  Maximization  of  the  vertical  strain  (o!*)22  with  constraint  on  vertical  Young’s  modulus 

^2  ^  >  5  GPa.  Horizontal  and  vertical  symmetry  is  specified. 

The  resulting  topologies  are  shown  in  Figs.  9,  10  and  11,  and  their  effective  properties  are 
shown  in  Table  2. 

To  overcome  the  positive  thermal  expansion  of  other  surrounding  materials,  we  seek  to 
maximize  the  contraction  force,  i.e.,  minimize  the  isotropic  thermal  stress  coefficient  as  in 
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Figure  9:  Examples  (f):  Optimal  microstructure  for  minimization  of  the  isotropic  thermal 
stress  coefficient  The  white  regions  denote  void  (phase  0),  the  filled  regions  consist  of 
Invar  (phase  1)  and  the  cross-hatched  regions  consist  of  Nickel  (phase  2). 


Example 

Objective 

cWjcW 

GPa 

kPa/K 

Invar 

0.8 

150 

0.31 

175 

1/0 

Nickel 

(Max. 

13.4 

200 

0.31 

3890 

0/1 

15 

0.06 

-78 

0.60/0.28 

(g)  Fig.  10 

10.0/-1.6 

9.2/8.8 

-0.75/-0.80 

258/-210 

0.49/0.38 

(h)  Fig.  10 

5.4/-4.7 

67/30 

0.06/0.02 

392/-128 

0.60/0.30 

(i)  Fig.  11 

Max.  022^ 

23/35 

1. 1/5.0 

-.14/-.62 

2/174 

0.38/0.46 

Table  2:  Thermoelastic  parameters  for  optimal  microstructures  made  of  Invar  (phase  1)  and 
Nickel  (phase  2). 
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Figure  10:  Examples  (g)  (top)  and  (h)  (bottom):  Optimal  microstructures  for  minimization 
of  thermal  stress  coefficient  (top)  and  minimization  of  vertical  contraction  stress 
The  white  regions  denote  void  (phase  0),  the  filled  regions  consist  of  Invar  (phase  1)  and  the 
cross-hatched  regions  consist  of  Nickel  (phase  2). 


Figure  11;  Example  (i):  Optimal  microstructure  for  maximization  of  thermal  strain  in  the 
vertical  direction  022^.  The  white  regions  denote  void  (phase  0),  the  filled  regions  consist  of 
Invar  (phase  1)  and  the  cross-hatched  regions  consist  of  Nickel  (phase  2). 

example  (f).  The  obtained  isotropic  contraction  stress  of  example  (f)  is  =  — 78GPa. 
By  relaxing  the  isotropy  requirement  and  allowing  orthotropic  materials  the  directional  con¬ 
traction  stress  can  be  increased.  In  example  (g)  we  minimize  the  value  of  1^22  and  get  the 
effective  value  ^22  ~  —  128GPa.  Minimizing  the  value  of  gives  us  a  composite  which  for 
fixed  boundaries  has  high  contraction  force  (remember  that  the  thermal  stress  coefficient 
is  the  stress  in  a  material  constrained  at  the  boundaries).  If  we  want  to  maximize  the  contrac¬ 
tion  force  for  a  material  with  free  boundaries,  we  should  minimize  the  product  (£^*)2(q:*)22  as 
done  in  example  (h).  The  “free  boundary”  stress  of  example  (g)  is  (£',)2(q:*)22  =  — 14GPa, 
whereas  the  “free  boundary”  stress  of  example  (g)  is  (E,)2(Q!.)22  =  — 141GPa. 

If  we  want  to  maximize  the  expansion  stress  of  the  composite,  the  best  choice  would  be 
to  take  solid  nickel  material  both  for  the  isotropic  and  the  directional  cases. 

The  isotropic  negative  thermal  expansion  materials  in  examples  (a),  (b)  and  (f)  all 
have  positive  Poisson’s  ratios  (0.04,  0.52  and  0.06,  respectively),  showing  that  there  is  no 
mechanistic  relationship  between  negative  thermal  expansion  and  negative  Poisson’s  ratio. 

In  example  (i)  we  see  again  that  allowing  orthotropy  can  lead  to  high  directional  expan¬ 
sion  coefficients.  The  vertical  coefficient  (q:*)22  of  example  (i)  is  2.6  times  higher  than  for 
solid  nickel,  but  at  the  cost  of  a  low  vertical  Young’s  modulus  (2.5%  of  solid  nickel). 


5  Conclusions 


We  have  proposed  a  method  to  design  material  microstructures  with  extreme  thermoelastic 
properties.  The  optimization  procedure  has  been  shown  to  be  very  accurate  in  producing 
the  optimal  microstructures.  Indeed,  the  results  of  this  study  motivated  Gibiansky  and 
Torquato  to  improve  upon  the  29'year  old  Schapery-Rosen-Hashin  bounds  on  the  thermal 
expansion  of  three-phase  media.  Our  obtained  values  are  close  to  the  Gibiansky- Torquato 
bounds.  We  have  shown  that  extreme  thermal  expansion  behavior  can  be  obtained  but  at 
the  cost  of  a  low  bulk  modulus.  Therefore,  there  is  a  tradeoff  between  extremizing  thermal 
strain  coefficients  on  the  one  hand  and  ending  up  with  a  stiff  material  on  the  other.  We 
have  also  shown  that  extreme  directional  thermal  expansion  can  be  obtained  by  allowing 
anisotropy  of  the  composites. 

For  the  topology  optimization  method  in  general,  the  results  in  this  paper  shows,  that 
the  method  produces  designs  which  are  optimal  indeed. 

In  practice,  how  can  our  optimally  designed  materials  be  manufactured?  They  may  be 
fabricated  (with  cell  sizes  down  to  a  few  millimeters)  using  stereolithography  techniques  [20] 
or  using  surface  micromachining  techniques  (with  cell  sizes  down  to  50  micros)  as  seen  for 
materials  with  negative  Poisson’s  ratios  in  Larsen,  Sigmund  and  Bouwstra  [26].  Furthermore, 
it  will  be  interesting  to  examine  whether  the  lessons  learned  from  this  continuum  analyses 
can  be  exploited  to  optimally  design  and  synthesis  of  materials  at  the  molecular  level  (e.g. 
Baughman  and  Galvao  [3]). 

Finally,  we  note  that  the  method  is  applicable  to  design  of  smart  materials  (piezoelectric 
or  shape- memory-alloy  inclusions).  In  a  future  paper,  the  procedure  described  here  will 
be  used  to  find  the  structures  that  optimize  the  piezoelectric  properties  of  the  material  for 
use  as  actuators  or  sensors.  The  method  can  also  be  modified  to  handle  three-dimensional 
microstructures.  The  extension  to  three  dimensions  is  straight  forward,  but  computer  time 
will  increase  dramatically.  Extensions  to  three  dimensions  for  two  material  phases  have  been 
done  in  Sigmund  [43]  and  Sigmund  and  Torquato  [44]. 
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A  Homogenization  theory 


This  appendix  summarizes  the  important  equations  for  computation  of  the  effective  ther¬ 
moelastic  properties  of  a  periodic  inhomogeneous  material  using  the  homogenization  theory 
as  developed  in  Bensoussan  et  al.  [5]  and  Sanchez-Palencia  [37]. 

We  want  to  find  the  effective  homogenized  thermoelastic  tensors  C/*],,  and 
of  a  periodic  material  microstructure  described  by  the  rectangular  base  cell  Y.  Using  the 
symmetries  of  the  thermoelastic  tensors,  Cijki  =  Cjiki  =  Cijik  =  Ckiij,  Aj  =  Pji  and  =  aji, 
the  anisotropic  elasticity  tensor  has  nine  and  the  thermal  expansion  tensors  have  three 
independent  parameters. 

Only  considering  first  order  terms  in  the  asymptotic  expansion,  it  can  be  shown  that 
relations  for  obtaining  the  effective  elasticity,  thermal  stress  and  thermal  strain  tensors  can 
be  written  in  energy  forms  as 

=  7Xcp,«(o„-e?,(r))(4P-£W(x'’))<iy,  (i6) 

4*’  = 

where  the  displacements  fields  and  F  are  solutions  to  the  following  cell  problems;  find 

€  V  and  F  €  U,  such  that 


I 

Y  dyg  dy, 


=  {u  :  V  is  Y  —  periodic] , 


Y  dyi  dyj 


=  /  vt,  e  V, 

Jy  dyj 

=  {u  :  V  is  Y  —  periodic}. 


The  fluctuation  strains  and  e^^{V)  are  defined  through  the  strain-displacement 

relations  eW(x«)  =  lidx'i'ISy,  +  dxfldy,)  and  e^,(r)  =  \(dT,ldy,  +  dVJdy,)  and  eJW 
are  three  linearly  independent  test  strain  fields.  In  Eqs.  (16)  and  (17),  Cijki,  I3ij  and 
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are  the  locally  (dependent  on  y)  varying  stiffness,  thermal  stress  and  thermal  strain  tensors, 
respectively. 

As  test  strain  fields  we  choose  the  three  unit  tensors  (in  two  dimensions)  = 
(1,0, 0,0),  =  (0,1, 0,0)  and  =  (0,0, 1,0)  and  can  be  ignored  for  symmetry 

reasons. 

The  equilibrium  equations,  Eqs.  (17),  are  solved  using  the  finite-element  method.  The 
base  cell  is  discretized  by  finite  elements  and  solving  Eqs.  (17)  means  solving  a  finite-element 
problem  with  periodic  boundary  conditions  for  the  three  different  prestrain  cases:  horizontal 
unit  strain,  vertical  unit  strain,  shear  unit  strain  as  given  by  the  tensors  and  for  a 
thermal  strain  field  resulting  from  the  locally  varying  thermal  stress  tensor  /S,j. 

Having  discretized  the  base  cell  by  N  finite  elements,  the  integrals  for  the  evaluation 
of  the  homogenized  properties  in  Eqs.  (16)  can  be  evaluated  on  the  element  level  and  the 
effective  properties  can  be  written  as  the  sums 

^  e=l 

(18) 

where  T®  is  the  area  of  element  e. 

For  a  more  thorough  description  of  the  numerical  homogenization  procedure  and  finite- 
element  discretization,  the  reader  is  referred  to  the  numerical  works  of  Bourgat  [6]  and 
Guedes  and  Kikuchi  [14]. 

To  solve  the  linear  system  of  equations  we  use  a  so-called  Element- By-Element  Precon- 
ditioned-Conjugate-Gradient  solver  (EBE-PCG).  The  PCG  solver  is  described  in  Numerical 
Recipes  [33]  and  its  application  to  finite  element  problems  is  discussed  in  Hollister  and 
Riemer  [19],  who  use  the  method  for  the  microstructural  analysis  of  human  bone  structure 
discretized  by  up  to  one  million  finite  elements.  The  advantages  gained  by  using  the  EBE- 
PCG  solver  for  the  present  design  method  are  multiple.  The  EBE-PCG  solver  is  an  iterative 
solver  using  a  starting  guess  for  the  displacement  vector.  By  using  the  displacement  vector 
from  the  preceding  design  step,  computational  time  can  be  saved.  Furthermore,  we  do  not 
need  an  exact  solution  to  the  finite-element  problem  in  the  beginning  of  the  design  sequence, 
thus  we  can  stop  the  solver  when  an  approximate  solution  has  been  reached.  Only  in  the  final 
iterations  of  the  optimization  procedure  we  need  an  exact  solution  and  then,  the  convergence 
requirements  of  the  solver  can  be  made  stricter.  Another  advantage  of  the  EBE-PCG  solver 
is,  that  is  does  not  require  an  assembly  of  the  global  stiffness  matrix  thereby  saving  storage 
space,  in  fact,  it  is  only  necessary  to  store  two  element  stiffness  matrices,  implying  that  the 
stiffness  matrix  for  a  particular  element  can  be  calculated  using  Eq.  (4).  Finally,  the  EBE- 
PCG  solver  eliminates  problems  with  increase  of  bandwidth  due  to  the  periodic  boundary 
conditions.  Opposing  nodes  of  the  base  cells  are  simply  given  the  same  node  numbers  causing 
no  increase  in  computational  complexity,  due  to  the  element-by-element  nature  of  the  EBE- 
PCG  solver. 


32 


B  Sensitivity  analysis 

This  appendix  derives  all  the  sensitivities  which  are  necessary  to  solve  the  sequential  linear 
programming  problem  Eq.  (10).  The  sensitivities  can  be  found  directly  from  the  strain  fields 
already  computed  by  the  homogenization  procedure. 

The  sensitivities  of  the  effective  stiffness  tensor  in  Eq.  (18)  with  respect  to  design  variables 
x\  and  xl  of  element  e  can  be  found  as 

^  X,  -  e<;Xx‘^))dY\  m  =  1, 2.  (19) 


where  it  was  used  that  the  test  fields  are  independent  of  the  design  variables  (i.e. 

def''>ldxl  =  =  0). 

The  sensitivities  of  the  local  stiffness  tensor  are 

^=(l-i|)CW  and  ^  =  xf(-C.<;2, +  C«).  (20) 

Sensitivity  of  the  effective  thermal  stress  tensor  with  respect  to  design  variable  xl  in 
Eq.  (18)  can  be  found  as 

^  =  “?X.  “  44r))(4i"’  -  (21) 

where  it  was  used  that  the  thermal  test  field  a,-;  is  independent  of  design  variable  x^  (i.e. 
daijfdxl  =  0). 

As  aij  is  dependent  on  design  variable  xl  (i.e.  da^jjdxl  —  sensitivity  of 

the  effective  thermal  tensor  with  respect  to  design  variable  xl  in  Eq.  (18)  has  an  extra  term 


dxl  Y  Jy‘  dx 


Sensitivities  of  the  effective  thermal  strain  tensor  aj*  with  respect  to  design  variables  a:® 
and  xl  are  found  by  differentiating  Eq.  (16)  as  follows: 

^  m  =  1,2.  (23) 


Finally,  the  sensitivities  of  the  volume  fractions  Eqs.  (5)  with  respect  to  design  variables 
x\  and  xl  are 

dc^^^/dxl  =  (1  -  xDY^jY,  dc^^^ldxl  =  x\Y^IY, 

dc^^^ldxl  =  x\Y^lY.  (24) 
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Figure  12:  Idealization  of  the  cross  microstructure  from  Fig.  8. 


C  Beam  model  of  thermal  actuator  example 

To  check  the  results  obtained  with  the  topology  optimization  algorithm,  this  appendix  calcu¬ 
lates  the  effective  thermal  strain  tensor  for  the  “maximization  of  vertical  expansion”  (Fig.  8) 
example  (e)  using  a  simple  beam  model. 

The  cross-like  topology  from  Fig.  8  can  be  idealized  as  seen  in  Fig.  12.  In  the  interpre¬ 
tation  we  assume  that  each  of  the  four  legs  are  slender  beams  of  height  /i,  width  ru  =  1 
and  length  I  —  \/5/2.  The  four  beam  legs  are  joined  together  at  the  center  of  the  cell  and 
connected  to  neighboring  cells  with  ideal  moment-free  hinges. 

The  thermal  strain  coefficient  of  the  beam  model  can  be  found  by  calculating  the  bending 
and  elongation  of  each  “leg”  and  projecting  these  displacements  to  the  horizontal  and  vertical 
directions.  Using  Bernoulli- Euler  beam  theory,  the  strain  across  the  thickness  of  a  beam 
varies  with  the  distance  from  the  neutral  axis  as  =  y/r,  where  is  the  strain  in  the 
x-direction  (horizontal)  and  r  is  the  radius  of  curvature.  Using  Hooke’s  law,  the  stress  in  the 
^-direction  of  the  beam  is  <7^  =  Eex  -  Eyjr.  The  thermal  expansion  of  the  two  layers  can 
be  seen  as  an  applied  thermal  stress  to  the  two  layers,  i.e.,  o-jp)  =  Edi  {i  =  1,2).  To  find  the 
beam  moment-curvature  relationship,  resulting  moments  about  the  transverse  beam-axis  are 
summed  to  give 

fh/2  rh/2 

w  /  (JTydy  =  It!  /  a^ydy  (25) 

J-h/2  J-h/2 


/  \  /A, 

/  Eoi^^’ydy  —  w  Ea^'^^ydy  =  w  E—dy. 
J-hj2  7o  J~hi2  r 


Eq.  (26)  can  be  evaluated  as 


w\Eh?a^^^y  —  w\Eh^a^^\  =  wE-^. 
8  8  12r 


(26) 


(27) 


Using  the  approximation  d^yfdx^  ~  I/?",  the  differential  equation  for  bending  of  the  beam 
can  be  found  as 

<^^2/  _  1  _  3(a<^) -0:^2)) 


dx'^ 


2h 
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Integrating  Eq.  (28)  twice  with  respect  to  x,  we  get  the  vertical  displacement  at  the  tip  of 
the  cantilever 


Uy 


Ah 


(29) 


The  displacement  in  the  ar-direction  is  simply  the  average  elongation 


(30) 


Rotating  the  leg  45  degrees,  the  resulting  deflection  of  one  leg  in  the  horizontal  and 
vertical  directions  are 


Uh  =  +  Uy 


and  =  ~{-u^  +  Uy). 


(31) 


The  effective  thermal  strain  coefficients  are  the  relative  displacements  of  the  whole  cross 
in  the  horizontal  and  vertical  directions,  respectively 


ly  ly  h  \  Ah  J 


(32) 


Finally,  by  inserting  /  =  >/2/2,  we  get 


«!t'  =  ^(c">  +  o‘«)±|f(aW-aW). 


(33) 


The  thickness  of  the  legs  in  Fig.  8  is  estimated  to  h  —  0.18.  Inserting  this  value  together 
with  material  data  in  Eq.  (33),  we  get  =  —21  and  022^  =  32. 

Again  using  simple  beam  analysis,  the  effective  bulk  modulus  can  be  found  as  Ri 
IQEwh^.  Inserting  ft  =  0.18  and  in  =  1,  we  get  0.085. 


35 


The  DCAMM  reports  are  issued  for  early  dissemination  of  research  results  from  the  Department  of 
Energy  Engineering,  Department  of  Hydrodynamics  and  Water  Resources,  Department  of  Solid  Mechanics, 
Department  of  Structural  Engineering  and  Materials,  Department  of  Naval  Architecture  and  Offshore 
Engineering,  Department  of  Mathematics,  and  Department  of  Mathematical  Modelling  at  the  Technical 
University  of  Denmark.  These  technical  reports  are  as  a  rule  submitted  to  international  scientific  journals 
and  should  not  be  widely  distributed  until  after  the  data  of  publication.  Whenever  possible  reference 
should  be  given  to  the  final  publications  and  not  to  the  DCAMM  reports. 


This  report  or  any  part  thereof  may  not  be  reproduced  without  written  permission  of  the  Danish  Center 
for  Applied  Mathematics  and  Mechanics. 


